options(rgl.useNULL=TRUE)

require(shinyRGL)
require(rgl)

#---------------------------Total Energy Consumption----------------------------#
energia<-read.csv('data/consumo_energia.csv',header=TRUE,sep=';')

totec<-energia[,-3];vnames<-names(totec[,3:12])
#totec<-reshape(totec,v.names=vnames,timevar='year',idvar='Divipola',direction='wide')


#################################################################################
#--------------------Per Capita Total Energy Consumption------------------------#

pob<-read.csv('data/poblacion.csv',header=TRUE,sep=';')

pob<-reshape(pob,varying=list(5:40),sep='',idvar='Divipola',v.names='pob',
             timevar='year',times=1985:2020,direction='long')

# Sort the data by the municipal numeric identifier 'divipola'
totec<-totec[ order(totec$Divipola,totec$year),]
pob<-pob[ order(pob$Divipola,pob$year),]
#identical(pob$Divipola,totec$Divipola)

year<-c(1996:1997,1999:2011)
pob<-pob[pob$year>=1996 & pob$year<=2012,]
pc<-totec


pc[,vnames]<-totec[,vnames]/pob[,'pob']


pc[is.na(pc)==1]<-0
for(i in 1:ncol(pc)){
    pc[pc[,i]=='Inf',i]<-0    
}

# Drop unconnected regions
totec<-subset(totec, !(totec$Divipola %in% c(5475,5873)))
pc<-subset(pc, !(pc$Divipola %in% c(5475,5873)))

#################################################################################
#-----------------------------Antioquia's Shape---------------------------------#

library(spdep)
library(maptools)

# Load Antioquia's shape file
ant<-readShapePoly('data/antioquia.shp')


# Set rownames, namely, change feature identifier in spatial object
row.names(ant)<-as.character(ant$ID_ESPAC_2)

# Check rownames, to confirm that the preceeding had the desired effect
#identical(row.names(ant),sapply(ant@polygons, function(x) slot(x,"ID")))

# Re-order
ant<-ant[order(row.names(ant)),]

# Indeed, the reordering worked
#row.names(ant)
ant<-subset(ant, !(ant$ID_ESPAC_2 %in% c(5475,5873)))

#################################################################################
#-----------------------------Travel Times Matrix-------------------------------#

odh<-read.table('data/od-matrix-h.csv',sep=';')

#identical(pc$Divipola,totec$Divipola,as.numeric(ant$ID_ESPAC_2),as.numeric(rownames(odh)))

#################################################################################
#------------------------------Market Potential---------------------------------#

#W<-as.matrix(odh);colnames(W)<-rownames(W)

#################################################################################
library(ks)

library(shiny)

# Define server logic for random distribution application
shinyServer(function(input, output, session) {
    
    data <- reactive({
        db<-pc
        
        x0 <- rowSums(cbind(0,db[db$year==input$year[1],input$var]))
        y0 <- rowSums(cbind(0,db[db$year==input$year[2],input$var]))
        
        x <- x0[x0>0 & y0>0]
        y <- y0[x0>0 & y0>0]
        
        rm(x0);rm(y0)
        
        if(input$log==TRUE){
            x <- log(x)
            y <- log(y)
        }
        
        if(input$std=='mean'){
            x <- x/mean(x)
            y <- y/mean(y)
        } else if(input$std=='median'){
            x <- x/median(x)
            y <- y/median(y)
        }
        
        list(x=x,y=y)
    })
    
    data1 <- reactive({
        db <- pc
        
        y <- db$year
        x <- rowSums(cbind(0,db[,input$var1]))
        
        y <- y[x>0]
        x <- x[x>0]
                
        if(input$log1==TRUE){
            x <- log(x)            
        }
        
        if(input$std1=='mean'){
            x <- x/mean(x)        
        } else if(input$std1=='median'){
            x <- x/median(x)        
        }
        
        list(x=x,year=y)
    })
    
    fhat <- reactive({
        
        H.pi<-Hpi(cbind(data()$y,data()$x), binned=TRUE)*3 ## a is a matrix of x,y,z points  
        fhat<-kde(cbind(data()$y,data()$x), H=H.pi,gridsize=100)
        return(fhat)
    })
        
    output$plot1 <- renderPlot({
        
        w <- density(data()$x,na.rm=TRUE,bw='SJ')
        z <- density(data()$y,na.rm=TRUE,bw='SJ')
        
        xlim <- c(min(w$x,z$x),max(w$x,z$x))
        ylim <- c(min(w$y,z$y),max(w$y,z$y))
        
        
        plot(w, main=paste(input$var,collapse='-'),ylab="Densidad",xlim=xlim,ylim=ylim)
        lines(z,col='blue')
        legend('topright',legend=c(input$year[1],input$year[2]),col=c('black','blue'),lty=c(1,1))
        
    })
    
    output$plot2 <- renderWebGL({
        
	#plot(fhat, display="persp", cont=seq(5,95,by=10),main=paste(input$var,collapse='-'),
        #     ylab=input$year[2],xlab=input$year[1],border=1,thin=2)
        
        rgl.viewpoint(theta=0,phi=-70,fov=50)
        persp3d(fhat()$eval.points[[1]],fhat()$eval.points[[2]],fhat()$estimate,
                col='green',front='line',back='line',alpha=0.8,fog=TRUE,axes=FALSE,
                xlab=input$year[1],ylab=input$year[2],zlab='Density',lwd=1)
        axes3d(c('x--','y--','z'),tick=FALSE,labels=FALSE)
    })
    
    output$plot3 <- renderPlot({

        plot(fhat(), display="filled.contour2", cont=seq(5,95,by=10),main=paste(input$var,collapse='-'),
             ylab=input$year[2],xlab=input$year[1])
        plot(fhat(), display="slice", cont=seq(5,95,by=10),drawlabels=TRUE,
             labcex=0.6,add=TRUE)        
        abline(0,1)
    })
    
    output$plot4 <- renderPlot({
        #data<-get(input$data)
        
        z<-tapply(data1()$x,data1()$year,sd)
        
        par(mfrow=c(2:1))
        boxplot(data1()$x~data1()$year,ylab=paste(input$var,collapse='-'),main='Boxplot')
        plot(z~names(z),t='b',ylab=paste(input$var,collapse='-'),xlab='',main='Standard Deviation')
    })


})
